Stochastic Simulation of a finite-temperature one-dimensional Bose-Gas: from 

Bogoliubov to Tonks-Girardeau regime 
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We present an ab initio stochastic method for calculating thermal properties of a trapped, ID 
Bose-gas covering the whole range from weak to strong interactions. Discretization of the problem 
results in a Bose-Hubbard-like Hamiltonian, whose imaginary time evolution is made computation- 
ally accessible by stochastic factorization of the kinetic energy. To achieve convergence for low 
enough temperatures such that quantum fluctuations are essential, the stochastic factorization is 
generalized to blocks, and ideas from density-matrix renormalization are employed. We compare 
our numerical results for density and first-order correlations with analytic predictions. 
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For a long time the homogeneous ID gas was a rather 
academic toy model, mainly of interest for mathematical 
physics since it is analytically tractable using the Bethe 
Ansatz Q, 0, Q • Recent successes in trapping and cool- 
ing bosonic atoms have lead however to a number of ex- 
periments culminatingin the observations of Tonks-gas 
properties in lattices Q and harmonic traps [j| • These 
developments have also lead to a rapidly growing theoret- 
ical interest in trapped ID gases, which can be described 
by the grand canonical Hamiltonian 



K 



dx ffi (x) 



ijj(x) 



dx^>' (x)^' {x)il){x)ip(x). (1) 



Here m is the mass of the bosons, <?id — the ID interac- 
tion constant, V{x) — the (harmonic) confining potential, 
and fi — the chemical potential. 

Physics of degenerate quantum gases in one dimension 
is dominated by quantum fluctuations and thus differs 
considerably from that in higher dimensions. Homoge- 
neous bosons do not form a condensate in ID 0. At 
any T > first-order correlations decay exponentially 
with distance; a power-law decay occurs only at T = 
0,0. The most peculiar feature of ID interacting bosons 
is that they become less ideal with decreasing density, 
n: the bosons start avoiding each other as if they were 
noninteracting fermions. This regime, called the Tonks- 
Girardeau gas is reached if the Tonks parameter, 
7 = mgii)/h 2 n, is large compared to unity. 

Experimentally, a quasi- ID regime can be realized in 
highly anisotropic traps, where the radial motion is con- 
fined to zero point oscillations and thus effectively elimi- 
nated. In the presence of a ID trapping potential, exact 
solutions are no longer available. Local properties can be 
derived in a local density approximation from the Bcthc- 
Ansatz solutions ^3 ■ The finite size of the trapped sys- 
tem results in an energy gap, which allows for the for- 
mation of a quasi-condensate in the weakly interacting 
regime 7 < 1. In this limit the quasi-condensate den- 
sity can be obtained from the Gross-Pitaevskii equation, 



and correlations can be calculated within Bogoliubov and 
Thomas- Fermi approximations [T3. IhH IH| . 

Long-distance correlations in the homogeneous gas 
can be described for all values of 7 using an effective 
harmonic-fluid (Luttinger liquid — LL) model 0,0, Eli 
which is analytically tractable. At T = it yields an 
algebraic decay of phase correlations ~ l/x 1 / 2K : where 
K > 1 is the so-called Luttinger parameter. In the case 
of a trap the ID gas can also be related to a LL model 
[l?! and solved within a local density approximation . 
In a trap, the power law sets in at temperatures of the 
order of the trap frequency, k^T w hw. 

We here present a stochastic method to numerically 
calculate the quantum properties of a trapped finite- 
temperature ID Bose gas. We are particularly inter- 
ested in the transition from the Gross-Pitaevskii to the 
Tonks-Girardeau regime as manifested in the behavior 
of the density and first order correlations. So as to 
achieve convergence for temperatures as low as the char- 
acteristic phase temperature T p h introduced in |13| |. we 
employ an extension of the stochastic Gutzwiller ap- 
proach for the grand canonical density operator [ljj, [2(J • 
Like the stochastic Gutzwiller Ansatz, our method is 
based on a stochastic factorization of the kinetic energy 
of a discretized Hamiltonian, which is however applied 
to blocks of sites. Furthermore, the blocks are treated 
by methods inspired by density matrix renormalization 
group (DMRG) techniques jU • In the strong interaction 
regime the kinetic energy is effectively a perturbation, 
which makes our approach more suitable than, e.g., the 
stochastic Hartree approach 22], (which has been very 
successful in describing the occupation number statistics 
of a weakly interacting BEC) or phase-space techniques 
of quantum optics in imaginary time |23l 12 1| . Recent 
extensions of the latter involving stochastic gauges [2{| 
also allow access to the regime of strong interactions, yet 
seem to be limited to higher temperatures. 

We start from deriving a consistently discretized ver- 
sion of Q). Let equidistant grids in position- and quasi- 
momentum spaces be introduced, with grid constants Ax 
and Ap, respectively: AxAp — 2-kH/M, M being the 
number of grid points. This is equivalent to putting the 
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system into a box of size L = MAx with periodic bound- 
ary conditions and restricting the quasi-momentum to 
an interval of length P = MAp. The Bose field with 
discrctized modes dj corresponding to wave numbers 
kj = Pj/H is related to local bosonic operators via the 
discrete Fourier transformation: (j, I = 0, • • • ,M — 1) 

Af-l 

^TtfE^** 3 ' MJWii. (2) 
In terms of the discretised fields Hamiltonian (JJJ reads 

m 13 (3) 

i i 
For symmetric and sufficiently small Ax one finds 

which only contains nearest-neighbor hopping. Im- 
portantly, the kinetic energy matrix scales like Ax~ 2 , 
whereas the interaction — like Ax . 

The discretized grand canonical Hamiltonian is equiva- 
lent to a Bose- Hubbard model (BHM) with effective hop- 
ping J = ti 2 /2mAx 2 , effective on-site interaction U — 
gro/ 'Ax, and effective chemical potential /ibh = M — 2 J. 
The scaled hopping can be expressed in terms of the 
Tonks parameter at the trap center, J/U = 2/^n(0)Ax. 
For fixed parameters, the discretized model is identical to 
the continuous one for Ax — > 0. In this limit always the 
weak coupling case J/U — > of the corresponding Bose- 
Hubbard model is recovered, which prohibits an insulat- 
ing phase. Thus the ID gas corresponds to a superfluid 



phase of the BHM close to the line ^bh = —2 J. However, 
for finite values of Ax, which should be less than n(0) _1 
and small enough in order to stay outside the n = 1 in- 
sulator lobe, a strong interaction limit J/U <g; 1 is still 
possible for large values of 7. 

Our goal is to numerically evaluate the grand canoni- 
cal statistical operator p = Z^ 1 cxp(— (3K). To gain ac- 
cess to lower temperatures, we generalize the stochas- 
tic Gutzwiller Ansatz [2(j, applying it to blocks of sites. 
Namely, assume for purposes of explanation, that the to- 
tal grid of M sites is divided into two blocks of M/2 
sites each. Hamiltonian J3J may be written as, /C 
H1+H2+H12, where H12 = - 2m ' Aa .s ^V/2 & m/2+i + h.c, 

while Hi (Ti.2) acts in the Hilbert space of the first (sec- 
ond) block. Applying the Trotter decomposition and the 
split approximation, we have, 

R -m ra L-dTH le -dTH 2e -dTH 12 \ A ^ A^l, (5) 

where dr — (3/ A is the step in imaginary time. Let rj T 
be A independent complex Gaussian stochastic variables, 
one per each time step r, such that 777 = 0, rj 2 = 0, 

and \rj T \ 2 — h 2 dr /2m Ax 2 . (r here is regarded as an in- 
dex numbering steps rather than a continuous variable.) 
Then, 

e -drH 12 _ j _ dr? ^ i2 _ 

(1 - VtI>m/2 + Vrhl/2)(l - Vrb\f /2+1 + Vrhl/2+l), 

(6) 

and the evolution operator is "stochastically factorised" : 
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lib 



-drH 



2 i 1 - ^M/2+l + vfb M /2+l) 



(7) 



(Recall that r\ T are independent for different r.) General- 
isation to arbitrary number of blocks is straightforward. 

Introducing blocks we however face another problem: 
the blocks themselves should be sufficiently small to han- 
dle them numerically. Since the block Hilbert spaces 
are quite large even when restricting ourselves to small 
particle numbers, some truncation must be introduced. 
Choosing naively the lowest eigenstates of the block 
Hamiltonians is dangerous because they do not neces- 
sarily correspond to the lowest eigenstates of the whole 
system. A fix to this problem is well known from the 
DMRG techniques [2]]]. Following those ideas, we firstly 



build a thermal p-matrix of a block with one additional 
site at each end (environment). The most important 
states are then found by tracing out the environmental 
degrees of freedom, diagonalising the block p-matrix thus 
obtained, and finally taking the states corresponding to 
largest eigenvalues. 

Not unexpectedly, the convergence of the stochastic 
factorization method becomes worse when Ax or the tem- 
perature are decreased. Ax must be chosen small enough, 
for the discrete system to approximate well the continu- 
ous one; the number of particles per site must be small 
compared to unity 0, |2g . Reducing Ax leads to larger 
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quasi-momentum cut-off and thus to larger noise in the 
simulation of the kinetic energy. With the block factor- 
ization method and Aa; < l/n(0) we were able to reach 
a temperature of k^,T = fku, corresponding to T rs T p h- 
Going to even lower temperature is mainly a question of 
a suitable Hilbert-space truncation. Block artifacts, i.e. 
effects of badly chosen block states, although much re- 
duced by use of the environment, show up especially at 
low temperature. 

In figures ^ and |3 numerical results for the density and 
the first order correlations in a ID trap are shown for 
k^T = fku. The densities are compared to the predic- 
tions obtained from the Bethe-Ansatz solutions of Yang 
and Yang using a local density approximation [i.e. re- 
placing n by fi(x) = fj, — V(x)]. Also shown is the Tonks- 
gas limit, 7 = 00, at the given temperature, which is ob- 
tained using the mapping to a free Fermi gas. Apart from 
the case of 7 = 0.8 where block artifacts are still present, 
there is a very good agreement with the prediction of the 
Yang- Yang theory in local density approximation. The 
latter becomes invalid close to the edges of the gas, and 
thus larger deviations of the numerical simulation from 
the Yang- Yang theory occur. 

The Yang- Yang solution unfortunately giving no infor- 
mation about the correlations in the system, we compare 
the corresponding numerical results to different predic- 
tions valid either in the weak or strong interaction limits. 
The weak-interaction limit is described by a Bogoliubov 
approximation. Since the temperature is rather low we 
here have not taken into account thermal depletion of the 
quasi-condensate. In the opposite limit, 7 — > 00, correla- 
tions can be calculated by mapping impenetrable bosons 
hi, b\ to fermions Cj, cj via hi — Ylj<iO- ~ 2c] Cj) Cj, which 
leads to the expression for first-order correlations found 
by Efetov and Larkin |]| 

(SJS i )=Det(g«), i<j. (8) 

g y is a (j — i) x (j — i) matrix with elements (g^) n ,m = 
(ajj c m +i) — S nm +i/2, with n and m running from i to 
Fig. [3 compares the simulated correlations with the Bo- 
goliubov and Efetov-Larkin predictions. The expected 
transition from the Bogoliubov to the Efetov-Larkin be- 
havior in the Tonks limit with increasing 7 is clearly seen. 

To see the asymptotic behavior of the phase correla- 
tions more clearly, in fig. [21 we plot the first-order corre- 
lations normalized to the density (^(xy^f(—x))/n(x) for 
7 = 4 on a logarithmic scale. We show numerical results 
for three different temperatures, k^Tjhio = 1,2,3. For 
the lowest temperature we were able to reach in our sim- 
ulation, deviations from the pure exponential decay char- 
acteristic of higher temperatures can already be seen for 
intermediate distances. This is consistent with predic- 
tions of the Luttinger- liquid model 0] , namely, that the 
asymptotic behaviour of the correlations changes from 
exponential to a power law if the thermal energy k^T 
becomes much smaller than the trap energy fuv. (The 
spatial resolution of our simulations is insufficient to see 
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FIG. 1: Particle density of the Bose gas in a trap at k^,T = huj 
for different interaction strengths. Dots with error bars: 
stochastic simulation, solid line: prediction from Yang and 
Yang within the local density approximation, dashed line: 
Tonks (fermionization) limit. 



the short-distance behavior of the correlations which is 
not described by the LL model.) In order to go to 
even lower temperatures, we have also employed a finite- 
temperature variant of the DMRG method [27j . This 
method works well for the complementary regime of very 
low temperatures fcsT//kj < 0.1. The results show a 
smooth continuation of the behavior found in the present 
paper and will be discussed elsewhere. 

In summary we have presented a numerical method 
for the stochastic simulation of thermal properties of one 
dimensional Bose gases covering the whole range from 
the Bogoliubov regime of weak interactions to the Tonks- 
Girardeau regime of strong interactions. The method is 
based on a generalization of the stochastic Gutzwiller 
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Ansatz to blocks of sites in the discretized Hamiltonian. 
Combining the stochastic block factorization scheme with 
a suitable selection of the most important states of the 
block, we were able to reach the transition regime be- 
tween thermal and quantum dominated correlations. 
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FIG. 2: First-order correlations in the Bose gas for the same 
parameter regimes as in fig. Dots with error bars: stochas- 
tic simulation, solid line: Bogoliubov approximation, dashed 
line: Efetov-Larkin approximation. 
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